/* **************************************************** */
/* * Borwein Method for finding PI - Parallel Version * */
/* **************************************************** */

#include <gmp.h>
#include <pthread.h>
#include <semaphore.h>
#include <stdio.h>
#include <stdlib.h>

/* Number of iterations */
#define NUM_ITERS 2000
/* Size of array. For best results, NUM_ITERS % ARRAY_SIZE should be 0 */
#define ARRAY_SIZE 2000

/* Global array of mpf_t variables */
mpf_t a[ARRAY_SIZE+1], aSqrt[ARRAY_SIZE+1], b[ARRAY_SIZE+1], pn, pnp;
/* Global array of semaphores */
/* Note: First semaphore of each array is not used, but was maintained to facilitate indexing */
sem_t aSem1[ARRAY_SIZE+1], aSem2[ARRAY_SIZE+1], aSqrtSem[ARRAY_SIZE+1], bSem[ARRAY_SIZE+1];

/* Thread 1: for calculating 'a' coefficients */
void *aCalc(void) {
	int i;

	for(i = 0; i < ARRAY_SIZE; i++) {
		/* Calculate a[i+1] */
		mpf_ui_div(a[i+1], 1, aSqrt[i]);
		mpf_add(a[i+1], a[i+1], aSqrt[i]);
		mpf_div_ui(a[i+1], a[i+1], 2);

		/* Warn Threads 2 & 3 that a[i+1] is calculated */
		sem_post(&aSem1[i+1]);
		sem_post(&aSem2[i+1]);

		/* Calculate sqrt(a[i+1]) for use in next iteration */
		mpf_sqrt(aSqrt[i+1], a[i+1]);

		/* Warn Thread 2 that sqrt(a[i+1]) is calculated */
		sem_post(&aSqrtSem[i+1]);
	}

	pthread_exit(NULL);
}

/* Thread 2: for calculating 'b' coefficients */
void *bCalc(void) {
	int i;
	mpf_t bBottom;

	mpf_init(bBottom);

	/* First iteration uses values that are always ready when this thread starts, hence sem_wait is not used */
	mpf_add_ui(b[1], b[0], 1);
	mpf_mul(b[1], b[1], aSqrt[0]);
	mpf_add(bBottom, a[0], b[0]);
	mpf_div(b[1], b[1], bBottom);

	/* Warn Thread 3 that b[1] is calculated */
	sem_post(&bSem[1]);

	for(i = 1; i < ARRAY_SIZE; i++) {
		/* Wait for Thread 1 to finish a[i] */
		sem_wait(&aSem1[i]);

		/* Start calculating b[i+1] */
		mpf_add(bBottom, a[i], b[i]);
		mpf_add_ui(b[i+1], b[i], 1);

		/* Wait for Thread 1 to finish sqrt(a[i]) */
		sem_wait(&aSqrtSem[i]);

		/* Finish calculating b[i+1] */
		mpf_mul(b[i+1], b[i+1], aSqrt[i]);
		mpf_div(b[i+1], b[i+1], bBottom);

		/* Warn Thread 3 that b[i+1] is calculated */
		sem_post(&bSem[i+1]);
	}

	/* Free temporary variable */
	mpf_clear(bBottom);

	pthread_exit(NULL);
}

/* Thread 3: for calculating 'p' coefficients */
void *pCalc(void) {
	int i;
	mpf_t pBottom;

	mpf_init(pBottom);

	for(i = 0; i < ARRAY_SIZE; i++) {
		/* Wait for Thread 1 to finish a[i+1] */
		sem_wait(&aSem2[i+1]);

		/* Start calculating p[i+1] */
		mpf_add_ui(pnp, a[i+1], 1);
		mpf_mul(pnp, pnp, pn);

		/* Wait for Thread 2 to finish b[i+1] */
		sem_wait(&bSem[i+1]);

		/* Finish calculating p[i+1] */
		mpf_add_ui(pBottom, b[i+1], 1);
		mpf_mul(pnp, pnp, b[i+1]);
		mpf_div(pn, pnp, pBottom);
	}

	/* Free temporary variable */
	mpf_clear(pBottom);

	pthread_exit(NULL);
}

int main(int argc, char *argv[]) {
	int k, i, errorNumber = 0;
	pthread_t thread[3];

	/* Define precision */
	if(2 == argc)
		mpf_set_default_prec(atoi(argv[1]));
	else if(argc > 2) {
		printf("Usage: %s [prec]\n", argv[0]);
		exit(1);
	}

	/* Set a[0] = sqrt(2) */
	mpf_init_set_d(a[0], 2.0);
	mpf_sqrt(a[0], a[0]);
	/* Set sqrt(a[0]) */
	mpf_init_set(aSqrt[0], a[0]);
	mpf_sqrt(aSqrt[0], aSqrt[0]);
	/* Set b[0] = 0 */
	mpf_init(b[0]);
	/* Set p[i] = p[0] = 2 + sqrt(2) */
	mpf_init_set_d(pn, 2.0);
	mpf_add(pn, pn, a[0]);
	/* Initialise p[i+1] */
	mpf_init(pnp);

	/* Initialise the rest */
	/* Note that aSem1[0], aSem2[0], aSqrtSem[0] and bSem[0] are not used. */
	for(i = 1; i <= ARRAY_SIZE; i++) {
		mpf_init(a[i]);
		mpf_init(aSqrt[i]);
		mpf_init(b[i]);
		sem_init(&aSem1[i], 0, 0);
		sem_init(&aSem2[i], 0, 0);
		sem_init(&aSqrtSem[i], 0, 0);
		sem_init(&bSem[i], 0, 0);
	}

	/* Start process */
	for(k = 0; k < NUM_ITERS; k += ARRAY_SIZE) {
		/* Open Thread 1 for running ARRAY_SIZE times until global array is full */
		if((errorNumber = pthread_create(&thread[0], NULL, aCalc, NULL)) != 0) {
			fprintf(stderr, "%s: pthread_create failed with error code %d.\n", argv[0], errorNumber);
			exit(1);
		}
		/* Open Thread 2 for running ARRAY_SIZE times until global array is full */
		if((errorNumber = pthread_create(&thread[1], NULL, bCalc, NULL)) != 0) {
			fprintf(stderr, "%s: pthread_create failed with error code %d.\n", argv[0], errorNumber);
			exit(1);
		}
		/* Open Thread 3 for running ARRAY_SIZE times */
		if((errorNumber = pthread_create(&thread[2], NULL, pCalc, NULL)) != 0) {
			fprintf(stderr, "%s: pthread_create failed with error code %d.\n", argv[0], errorNumber);
			exit(1);
		}

		/* Wait for all threads to finish */
		for(i = 0; i < 3; i++) {
			if((errorNumber = pthread_join(thread[i], NULL)) != 0) {
				fprintf(stderr, "%s: pthread_join failed with error code %d.\n", argv[0], errorNumber);
				exit(1);
			}
		}

		/* Recycle array: Last elements become first and process is restarted */
		mpf_set(a[0], a[ARRAY_SIZE]);
		mpf_set(aSqrt[0], aSqrt[ARRAY_SIZE]);
		mpf_set(b[0], b[ARRAY_SIZE]);
	}

	/* Print final results */
	gmp_printf("Final a[n]: %.6Ff\n", a[0]);
	gmp_printf("Final b[n]: %.6Ff\n", b[0]);
	gmp_printf("Final p[n] (PI): %.6Ff\n", pn);

	/* Free everything */
	for(i = 0; i <= ARRAY_SIZE; i++) {
		mpf_clear(a[i]);
		mpf_clear(aSqrt[i]);
		mpf_clear(b[i]);
		if(i) {
			sem_destroy(&aSem1[i]);
			sem_destroy(&aSem2[i]);
			sem_destroy(&aSqrtSem[i]);
			sem_destroy(&bSem[i]);
		}
	}
	mpf_clear(pn);
	mpf_clear(pnp);

	return 0;
}